* Figure VII panels c)-d)
* 29 Oct 2021

clear all
set more off

grstyle init
grstyle set plain, nogrid horizontal
grstyle color background white

*Choose figure width:
local fig_w= 5
local res = `fig_w' * 600

import excel "$output_s/w4_test_of_model.xlsx", sheet("Sheet1") firstrow clear 

* Bar graph of cases - model vs prediction
gen occup_opt = ceil(case_opt_w4/2) //Quick way to separate cases
gen occup = ceil(case_real_w4/2)

label define occupation 1 "Mixing" 2 "Livestock Rearing Only" 3 "Wage Work Only" 4 "No Work"
label values occup* occupation

drop if occup == 4

preserve
	keep hhid occup*
	rename occup_opt occup2
	rename occup occup1
	reshape long occup, i(hhid) j(tag)
	label define tags 1 "Actual" 2 "Predicted"
	label values tag tags
	catplot tag occup, percent(tag) var1opts(label(labsize(small))) var2opts(label(labsize(small))) recast(bar) asyvars bargap(20) bar(1, color(pink*0.4)) bar(2, color(green*0.4)) ytitle("Percentage") name(bars, replace) yla(0(10)40, gmax)
	graph save bars "$output/Bars_cases_excl_no_work_w4.gph", replace
	graph export "$output/Bars_cases_excl_no_work_w4.png", as(png) replace
restore


set level 95 //Setting confidence interval


* Figure VII panel c)

	*Colored version
	twoway (lpolyci l_opt_w4 k_w4 if k_w4<16000, lcolor(green) fcolor(green%30) alcolor(%30) bw(800)) (lpolyci l_w4 k_w4 if k_w4<16000, lcolor(pink) fcolor(pink%30) alcolor(%30) bw(800)), xtitle("Capital at Wave 4 (BDT)", margin(t+2)) ytitle("Hours in Livestock Rearing") yla(0(200)800) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_l, replace)
	graph save Lpoly_l "$output/Lpoly_l_kw4.gph", replace
	graph export "$output/Lpoly_l_kw4.png", as(png) replace

	*TIFF Grayscale version:
	twoway (lpolyci l_opt_w4 k_w4 if k_w4<16000, lcolor(gs10) lp(dash) fcolor(gs10%40) alcolor(%40) bw(800)) (lpolyci l_w4 k_w4 if k_w4<16000, lcolor(gs4) fcolor(gs4%40) alcolor(%40) bw(800)), xtitle("Capital at Wave 4 (BDT)", margin(t+2)) ytitle("Hours in Livestock Rearing") yla(0(200)800) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_l, replace)
	graph export "$output/TIFF/Lpoly_l_kw4.tif", as(tif) replace width(`res')


* Figure VII panel d)

	*Colored version:
	twoway (lpolyci h_opt_w4 k_w4 if k_w4<16000, lcolor(green) fcolor(green%30) alcolor(%30) bw(800)) (lpolyci h_w4 k_w4 if k_w4<16000, lcolor(pink) fcolor(pink%30) alcolor(%30) bw(800)), xtitle("Capital at Wave 4 (BDT)", margin(t+2)) ytitle("Hours in Wage Work") yla(0(200)1400) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_h, replace)
graph save Lpoly_h "$output/Lpoly_h_kw4.gph", replace
graph export "$output/Lpoly_h_kw4.png", as(png) replace

	*TIFF Grayscale version:
	twoway (lpolyci h_opt_w4 k_w4 if k_w4<16000, lcolor(gs10) lp(dash) fcolor(gs10%40) alcolor(%40) bw(800)) (lpolyci h_w4 k_w4 if k_w4<16000, lcolor(gs4) fcolor(gs4%40) alcolor(%40) bw(800)), xtitle("Capital at Wave 4 (BDT)", margin(t+2)) ytitle("Hours in Wage Work") yla(0(200)1400) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_h, replace)
	graph export "$output/TIFF/Lpoly_h_kw4.tif", as(tif) replace width(`res')
